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Two  computationally  efficient  methods  for  superresolution  reconstruction  and  restoration  of  mi¬ 
croscanning  imaging  systems  are  presented.  Microscanning  creates  multiple  low-resolution  images 
with  slightly  different  sample-scene  phase  shifts.  The  digital  processing  methods  developed  here 
combine  the  low-resolution  images  to  produce  an  image  with  higher  pixel  resolution  (i.e.,  superreso¬ 
lution)  and  higher  fidelity.  The  methods  implement  reconstruction  to  increase  resolution  and  resto¬ 
ration  to  improve  fidelity  in  one-pass  convolution  with  a  small  kernel.  One  method  uses  a  small-kernel 
Wiener  filter  and  the  other  method  uses  a  parametric  cubic  convolution  filter.  Both  methods  are  based 
on  an  end-to-end,  continuous- discrete- continuous  microscanning  imaging  system  model.  Because 
the  filters  are  constrained  to  small  spatial  kernels  they  can  be  efficiently  applied  by  convolution  and 
are  amenable  to  adaptive  processing  and  to  parallel  processing.  Experimental  results  with  simulated 
imaging  and  with  real  microscanned  images  indicate  that  the  small-kernel  methods  efficiently  and 
effectively  increase  resolution  and  fidelity.  ©  2006  Optical  Society  of  America 
OCIS  codes:  100.6640,  100.3020,  100.2000. 


1.  Introduction 

Improvements  in  the  resolution  and  fidelity  of  digital 
imaging  systems  have  substantial  value  for  remote 
sensing,  military  surveillance,  and  other  applications, 
especially  those  for  which  a  large  field  of  view  is  desir¬ 
able  and  the  distance  to  objects  of  interest  cannot  be 
reduced.  Advances  in  optics  and  sensor  technologies 
offer  one  path  for  increasing  the  spatial  resolution  and 
fidelity  of  imaging  systems,  but  such  hardware  im¬ 
provements  can  be  costly  or  limited  by  physical  con¬ 
straints  or  both.1-2  Digital  image  processing  offers  an 
alternative  path  for  improving  image  quality. 

Digital  superresolution  reconstruction  and  restora¬ 
tion  methods  can  substantially  increase  resolution 
and  fidelity  of  an  imaging  system.  Reconstruction 
methods  increase  pixel  resolution  beyond  that  of  the 
physical  imaging  system  by  estimating  values  on  a 
finer  grid  or  lattice.  Restoration  methods  increase 
fidelity  beyond  that  of  the  physical  imaging  system  by 
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correcting  for  acquisition  artifacts  such  as  blurring, 
aliasing,  and  noise.  Superresolution  reconstruction 
and  restoration  can  be  combined  to  produce  images 
with  greater  resolution  and  higher  fidelity.3 

Typically,  superresolution  methods  produce  an  en¬ 
hanced  image  from  multiple  low-resolution  images 
with  either  different  sample-scene  phase  or  with  dif¬ 
ferent  blurring  functions.4-5  Microscanning  is  a  sys¬ 
tematic  approach  to  acquiring  images  with  slightly 
different  sample-scene  phases;  between  successive 
images  the  system  is  shifted  slightly,  either  in  a  pre¬ 
determined  pattern  or  in  a  random  pattern.  As  in 
most  superresolution  methods,  the  methods  devel¬ 
oped  in  this  paper  presume  global  sample-scene 
phase  shifts. 

Most  superresolution  methods  can  be  decomposed 
into  three  tasks:  registration,  reconstruction,  and  res¬ 
toration.  These  tasks  can  be  implemented  as  separate 
steps  or  in  one  or  two  combined  steps: 

Registration  is  the  process  of  orienting  several  im¬ 
ages  of  a  common  scene  with  reference  to  a  single 
coordinate  system.  Registering  images  to  subpixel  ac¬ 
curacy  is  a  crucial  factor  that  greatly  affects  superreso¬ 
lution  processing.  It  may  be  trivial  to  register  systems 
with  precisely  controlled  phase  shifting.  Popular  meth¬ 
ods  for  subpixel  image  registration  include  phase  or 
cross-correlation  registration6-7  and  gradient-based 
registration.8-9 
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Image  reconstruction  is  the  process  of  estimating 
image  value  at  arbitrary  locations  on  a  spatial  contin¬ 
uum  from  a  set  of  discrete  samples.  For  superresolu¬ 
tion  imaging,  the  aim  of  reconstruction  is  to  produce  an 
image  with  pixels  on  a  high-resolution  grid  with  uni¬ 
form  spacing  from  multiple  registered  images,  which 
may  have  pixels  that  are  irregularly  distributed. 
Popular  methods  for  reconstruction  include  nearest- 
neighbor  interpolation,  bilinear  interpolation,  cubic- 
spline  interpolation,  piecewise  cubic  convolution, 
and  cubic  optimal  maximal-order-minimal-support  (o- 
Moms)  interpolation.10  Kriging,  which  originated  in 
the  geostatistic  community,  is  a  popular  method  for 
interpolating  spatial  data.11-12 

Image  restoration  is  the  process  of  recovering  a 
more  accurate  image  of  a  scene  by  correcting  or  re¬ 
ducing  degradations  such  as  acquisition  blurring, 
aliasing,  and  noise.13  Basic  methods  for  image  resto¬ 
ration  include  deconvolution,  use  of  least-squares  fil¬ 
ters,  and  iterative  approaches.14 

Digital  processing  for  superresolution  is  an  area  of 
active  research.1'15  Approaches  based  on  nonuniform 
interpolation  are  the  most  intuitive.  Computational 
efficiency  is  their  major  advantage.  For  example,  Alam 
et  al.  proposed  a  superresolution  method  for  infrared 
images,  in  which  a  gradient-based  method  to  estimate 
shift,  a  weighted  nearest-neighbor  approach  to  align 
irregularly  displaced  low-resolution  images  to  form  a 
superresolution  image,  and  finally  a  Wiener  filter  for 
image  restoration  are  used.9  However,  nonuniform  in¬ 
terpolation  approaches  do  not  account  for  variations  in 
acquisition  conditions  for  the  low-resolution  images  or 
guarantee  the  optimality  of  the  end-to-end  imaging 
system.3  Tsai  and  Huang  investigated  superresolution 
imaging  in  the  frequency  domain,16  using  the  shift 
theorem  of  the  Fourier  transform  and  the  aliasing  re¬ 
lationship  between  low-resolution  images  and  the 
ideal  superresolution  image.  High-frequency  compo¬ 
nents  are  extracted  from  aliased  low-frequency  com¬ 
ponents.  Kim  et  al.  extended  this  approach  for  blurred 
and  noisy  images  and  developed  a  weighted  recursive 
least-squares  algorithm.17-18  Unfortunately,  Fourier- 
domain  methods  are  computationally  costly.  Super¬ 
resolution  methods  based  on  stochastic  theories 
employ  a  priori  knowledge  about  a  scene  and  noise. 
For  example,  Shultz  and  Stevenson19  proposed  a 
Bayesian  restoration  with  a  discontinuity-preserving 
prior-image  model,  and  Elad  and  Feuer5  and  Elad  and 
Hel-Or20  approached  superresolution  by  unifying  the 
maximum-likelihood,  maximum  a  posteriori,  and  pro¬ 
jection  onto  convex  set.  Stochastic  approaches  also 
may  have  high  computational  complexity  and  subop¬ 
timality  for  an  end-to-end  imaging  system.  Computa¬ 
tionally  efficient  spatial-domain  methods  have  been 
proposed.  Nguyen  and  Milanfar21  proposed  efficient 
block  circulant  preconditioners  for  solving  the  Tik¬ 
honov  regularized  superresolution  problems  by  using 
the  conjugate  gradient  method,  but  the  method  re¬ 
quires  perfect  shift  estimation.22  Farsiu  et  al.  proposed 
a  robust  method  that  defines  the  cost  function  of  reg¬ 


ularization  by  use  of  the  L1  norm  minimization.23  The 
method  is  robust  for  estimation  errors  of  shift  and 
noise,  but  the  system  model  accounts  only  for  acquisi¬ 
tion  in  the  end-to-end  system,  and  the  proposed  fast 
solution  is  iterative.  Even  though  interpolation  and 
restoration  are  done  simultaneously,  computational 
costs  for  real-time  applications  are  a  concern  with  it¬ 
erative  methods. 

In  this  paper  we  investigate  superresolution 
methods  for  microscanning  imaging  systems  that 
can  be  implemented  efficiently  with  small  convolu¬ 
tion  kernels.  The  approach  is  based  on  an  end-to- 
end  continuous- discrete- continuous  (CDC)  system 
model  that  accounts  for  the  fundamental  trade-off 
in  imaging  system  design  between  acquisition  blur¬ 
ring  related  to  optics,  detectors,  and  analog  circuits 
and  aliasing  owing  to  sampling.  The  system  model 
also  accounts  for  noise  associated  with  unpredict¬ 
able  signal  and  system  variations  and  quantization, 
for  misregistration  of  the  low-resolution  images, 
and  for  the  display  system.  The  approach  uses  one- 
pass  convolution  with  small  kernels  for  computa¬ 
tionally  efficient  reconstruction  and  restoration. 
This  approach  also  is  amenable  to  adaptive  process¬ 
ing  and  to  parallel  processing  with  appropriate 
hardware  support. 

The  rest  of  this  paper  is  organized  as  follows:  in 
Section  2  we  present  the  CDC  system  model,  introduce 
microscanning,  and  formulate  the  system’s  fidelity.  In 
Section  3  we  derive  superresolution  reconstruction  and 
restoration  filters,  including  the  optimal  Wiener  filter, 
a  spatially  constrained  Wiener  kernel,  and  a  paramet¬ 
ric  cubic  convolution  kernel.  In  Section  4  we  present 
experimental  results  for  images  from  a  simulation  and 
from  a  real  imaging  system.  In  Section  5  we  summa¬ 
rize  this  paper  and  describe  issues  for  future  work. 

2.  System  Model  and  Problem  Formulation 

In  this  section  we  present  the  CDC  system  model, 
describe  microscanning,  and  formulate  the  fidelity  of  a 
microscanning  imaging  system  based  on  the  system 
model.  The  system  model  is  introduced  in  the  spatial 
domain  of  the  image,  but  the  problem  and  fidelity  anal¬ 
ysis  are  formulated  in  the  Fourier  frequency  domain, 
so  spatial  convolution  can  be  considered  pointwise 
multiplication  of  transform  coefficients. 

A.  Continuous-Discrete-Continuous  System  Model 

The  superresolution  methods  developed  in  this  paper 
are  based  on  the  CDC  model  pictured  in  Fig.  1.  This 
imaging  system  model  is  relatively  simple,  yet  it 
captures  the  most  significant  degradations  in  typical 
imaging  systems:  linear  shift-invariant  blurring 
characterized  by  acquisition  point-spread  function 
(PSF)  h ;  aliasing,  which  is  due  to  sampling  a  contin¬ 
uous  function  on  a  uniform,  rectangular  lattice  I  I  I  ; 
additive  system  noise  e;  and  display,  characterized  by 
display  PSF  d. 

With  this  model,  a  single  low-resolution  digital  im¬ 
age  p  is  defined  mathematically  as 
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Fig.  1.  End-to-end  model  of  the  digital  imaging  process. 


p[m,  7l] 


s(x,  y)h(m  —  x,  n  —  y)dxdy 


pk[m,  ?z] 


s(x-xk,  y  -  yk ) 


+  e\m,  n], 


(1) 


X  hk(m  —  x,  n  — y)dxdy  +  ek[m,  ri\,  (3) 


where  [m,  n\  are  integer  pixel  indices  for  digital  im¬ 
age  p  and  (x,  y)  are  continuous  coordinates  for  scene 
s.  For  notational  convenience,  and  without  loss  of 
generality,  the  spatial  coordinates  are  normalized  in 
units  of  the  sampling  interval.  In  practice,  the  spatial 
extent  of  the  image  is  finite,  but  that  issue  is  not 
significant  for  the  following  analyses. 

Restoration  and  reconstruction  can  be  implemen¬ 
ted  in  one  step  by  convolution  of  image  p  with  filter 
PSF  f  to  produce  a  processed  image  of  arbitrary  reso¬ 
lution,  which  can  then  be  displayed.  Modeling  the  dis¬ 
play  as  a  linear  shift-invariant  process  characterized 
by  display  PSF  cL  yields  continuous  output  image  r. 


where  k  is  the  index  for  the  microscanning  image, 
(xk,  yk)  is  the  relative  shift,  lik  is  the  acquisition  PSF, 
and  ek  is  the  additive  system  noise. 

Image  acquisition  (with  blurring,  sampling,  and 
noise),  digital  processing  (for  registration,  reconstruc¬ 
tion,  and  restoration),  and  display  of  microscanned 
imaging  are  analyzed  more  easily  in  the  Fourier  fre¬ 
quency  domain,  regardless  of  whether  digital  image 
processing  is  performed  in  the  spatial  domain  or  in 
the  frequency  domain.  In  the  frequency  domain, 
the  Fourier  transform  of  microscanned  image  pk  [the 
transform  of  Eq.  (3)]  is 

pk(u,  o)=ns(B-|i,  v-  v)hk(u  -  p,  v  -  v) 

M-  v 

X  exp{— i2tt[(w  -  \x)xk  +  (u  -  v)yk]} 

+  ek(u,  v),  (4) 


X  d(x  -  x' ,  y -y')dx'dy'.  (2) 

B.  Microscanning 

Microscanning  is  the  process  of  generating  multiple 
images  from  a  common  scene  by  shifting  either  the 
scene  or  the  image-acquisition  system.  The  shifting 
can  be  performed  in  a  regular  pattern  or  an  irregular 
pattern.  Figure  2  illustrates  the  microscanning  pro¬ 
cess  for  a  sequence  of  images  pk,  k  =  0  .  .  .  K  -  1,  with 
an  unchanging  scene  shifted  between  images,  vari¬ 
able  blur  and  noise,  and  a  fixed  sampling  grid.  (Re¬ 
verse  shifting  of  the  sampling  grid  for  a  fixed  scene 
produces  the  same  images.)  Then  microscanned  im¬ 
age  pk  is 


Image 

Po 


Image 

Pk 


Fig.  2.  Microscanning  produces  multiple  images. 


where  a  circumflex  indicates  a  Fourier  transform.  In 
Eq.  (4)  the  frequency-domain  equivalent  for  spatial- 
domain  blurring  by  convolution  in  Eq.  (3)  is  pointwise 
multiplication  of  the  transform  coefficients  of  the 
scene  and  the  PSF.  The  frequency-domain  equivalent 
for  spatial-domain  sampling  is  the  double  sum,  which 
folds  the  transform  coefficients. 

The  microscanned  images  must  be  registered  rela¬ 
tive  to  one  another.  In  the  frequency  domain  the  reg¬ 
istered  and  combined  microscanned  images  are 

1  K-l 

p(u,  v)=  K  ^Pk(u,  o)exp{i2TT[M(xA  +  ak) 

+  v(yk  +  P*)]} 

1 

=  S(u  -  |X,  v-v) 

(jl  V 

K-l  yy 

x  2  hk(u  —  p,  v  —  v)exp[i2TT(wcq,  +  up*)] 

X  exp[t2iT(pXi  +  vy*)] 

1  K- 1 

+  ^  2  ek(u,  u)exp{i2TT[w(Xi  +  op) 

^  k=0 

+  v(yk  +  Pi)]},  (5) 
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where  (op,  (3/,)  is  the  registration  error  for  image  pk. 
Mathematically,  the  registered  images  are  combined 
by  addition.  If  registration  is  perfect,  then  op  =  (3* 
=  0,  and  each  microscanned  image  is  shifted  to  its 
proper  position  in  the  registered  image. 

In  the  Fourier  frequency  domain,  reconstruction 
and  restoration  of  the  scene  are  products  of  the 
Fourier  transform  of  registered  image  p  and  filter 
transfer  function  f.  Then  the  display  process  multi¬ 
plies  that  result  by  display  transfer  function  d: 


e2  =  J  J  [$*0,  v)~f(u>  v)d(u,  v) 

~  v)d*(u,  v)<hS:P(u,  v ) 

+  \f(u,  u)|2|d(zp  u)|2<Fp(m,  a)]dndu,  (9) 

where  d>p  is  the  power  spectrum  of  the  registered 
image  and  <T>V)  is  the  cross-power  spectrum  of  the 
scene  and  the  registered  image: 


r(u,  v)  =  p(u,  v)f(u,  v)d(u,  v ).  (6) 


C.  Fidelity  Analysis 

By  Rayleigh’s  theorem,  the  expected  mean  square 
error  (MSE)  of  the  CDC  imaging  system  for  an  en¬ 
semble  of  scenes  can  be  analyzed  in  either  the  spatial 
or  the  frequency  domain: 


e2  =  s 


| r(x,  y)  —  s(x,  y)  | 2  dxdy 


*(u,  v)-s(u,  u)|2dnduf.  (7) 


In  the  following  analysis  it  is  assumed  that  the 
power  spectra  of  the  scene  ensemble  and  the  noise  are 
known,  that  the  scene  and  the  noise  are  uncorrelated, 
that  coaliased  components  of  the  sampled  scene  are 
uncorrelated,  and  that  the  noise  between  images  is 
uncorrelated24: 


<bp(w,  v)  =  e  [\p(u,  u)|2] 


=  ^  S  S  (I\(b  -  p,  v  -  v) 


X 


v 

K-l 


E  hk(u  —  p,  v  —  v)e{exp[i2TT 


X  (wop  +  opb)]exp[i2iT 

X  (p**  +  vy*)]}  +  E  kh(u,  v), 
®s,p(u,  v )  =  e[s(u,  v)p*(u,  a)] 

-  1  K-l  A 

=  <bs(w,  aW  E  hk*(u,  u)e{exp[-i2TT 

k=0 

X  (wop  +  a|3*)]}.  (10) 


If  the  distribution  of  the  registration  errors  is 
known,  the  expressions  for  the  expected  MSEs  can  be 
analyzed  with  respect  to  those  errors.  For  example,  if 
op  and  |3,  are  uniformly  distributed  in  the  intervals 


1  1 

1  1  1 

2  Wp  2  Wx 

2 Wy’  2 Wy 

respectively,  then,  because 


e[s(w,  v)s*(u  —  p,  a  —  v)] 


d>s(u,  a)  (p,  v)  =  (0,  0) 
0  otherwise  ’ 


s[ej(u,  v)ek*(u,  a)] 


4*(“,  u)  j~k 

0  otherwise’ 


2Wy 

exp[i2-n(uak  +  v$h)~\WxWy&ak&$k 
i 

2  Wx  2  Wy 

=  sinc(w/Wt)sinc(a/W':v),  (11) 
the  components  of  the  expected  MSEs  are 


e[s(w  -  p,  a  -  v)ek*(u,  a)]  =  0,  (8) 

where  the  asterisk  denotes  complex  conjugation,  <3>s 
is  the  power  spectrum  of  the  scene,  and  is  the 
power  spectrum  of  the  noise.  For  convenience,  and 
without  loss  of  generality,  scenes  are  normalized 
such  that  the  mean  and  the  variance  are  0  and  1, 
respectively. 

The  expected  MSE  e2  for  the  CDC  imaging  system 
can  be  expressed  in  terms  of  the  scene  and  noise 
power  spectra,  the  acquisition  transfer  function,  the 
relative  shifts  and  registration  errors,  the  reconstruc¬ 
tion  and  restoration  filter,  and  the  display  transfer 
function: 


d>p(w,  a)  =  ^  sinc2(w/WI)sinc2(a/W’.y) 
x  E  E  chs(“  -  p,  v-v) 

X 


M-  V 

K-l 


E  hk(u  -  p,  a  -  v)exp[i2'iT(px<, 


+ 


+  vyh)\ 

d>sp(w,  a )  =  sinc(w/W*)sinc(a/Wy) 

A  1  K-l  „ 

X  d>s(w,  a)  K  E  hk*(u,  a).  (12) 


<!*,,  is  subject  to  relative  shifts  (xk,  yk).  For  instance,  if 
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the  number  of  microscanned  images  K  is  2,  <t>p  for 
relative  shifts  {(x0,  y0)  =  (0,  0),  {xx,  yx)  =  (0,  0.5))  is 
the  same  as  for  relative  shifts  {(x0,  y0)  =  (0.5, 
0),  (xj,  yx)  =  (0.5,  0.5)}  but  is  different  from  that  for 
relative  shifts  {(x0,  y0)  =  (0,  0),  (x1;  yj  =  (0.5,  0)}  and 
for  {(x0,y0)  =  (0,  0.5),  (xj,^)  =  (0.5,  0.5)}.  Filters  de¬ 
rived  with  respect  to  (f>;,  can  vary,  depending  on  the 
relative  shift  pattern  of  microscanned  images. 

Fidelity25  is  a  normalized  measure  of  image  quality 
based  on  the  MSE: 


3?  =  1 - .  (13) 

J  J  <Z>s(u,  u)dudu 

The  greatest  fidelity  possible  is  1  when  the  MSE  is  0. 
In  Section  3  we  derive  superresolution  reconstruction 
and  restoration  filters  that  maximize  fidelity  2F  (or 
equivalently  minimize  the  MSE,  e2). 

3.  Superresolution  Reconstruction  and  Restoration 

The  techniques  developed  in  this  section  are  designed 
for  performance  of  reconstruction  and  restoration. 
The  methods  take  as  input  multiple  images  that  have 
been  registered  and  produce  a  single  output  image 
with  high  fidelity  and  superresolution. 

A.  CDC  Wiener  Filter 

For  comparison  of  performance,  it  is  useful  to  derive 
the  optimal  CDC  Wiener  filter.24  Denoting  MSE  as  a 
functional  of  the  filter  transfer  function  f(u,  v )  yields 

e2(/)  =  J  J  L(u,  v,  f)dudv,  (14) 

where 

L(u,  v,  f)  =  <bs(w,  v )  - f(u ,  v)d(u,  v)d>S:P(u,  v ) 
-f\u,  v)d*(u,  o)<J \p*(u,  v ) 

+  | f(u,  u)|2|d(w,  U)|2d>p(w,  u).  (15) 

The  optimal  filter  must  satisfy 

=f*(u ,  v)\d(u,  u)|  %(u,  v)  —  d(u,  v)d>StP*(u,  v) 

dt 

=  0,  (16) 
so  the  optimal  filter  is 

ffy.pfy,  v )  d*(u,  v)  ®s,p(u,  v) 

fw{U,  V )  =  — -  7^ - f2  =  ^ - , 

T>p(w,  v)  |  d(u,  u)|  d>p(u,  v ) 

(17) 

where  rfyand  (1>S  p  are  introduced  to  incorporate  the 


effects  of  the  display  device  on  the  image  (and  are  used 
in  the  derivations  of  the  small  kernels  in  the  following 
sections).  The  CDC  Wiener  filter  cannot  be  imple¬ 
mented  practically  by  means  of  spatial  convolution  be¬ 
cause  it  is  continuous  and  its  support  is  the  full  extent 
of  the  image.  As  described  in  Subsection  3.B,  one  can 
reduce  the  computational  costs  of  superresolution  re¬ 
construction  and  restoration  by  constraining  the  spa¬ 
tial  support  of  the  filter  to  a  small  kernel. 

B.  Small-Kernel  Wiener  Filter 

The  derivation  of  small-kernel  Wiener  filter  fc  is  con¬ 
ditioned  on  constraints  imposed  on  its  spatial  sup¬ 
port.  The  support  of  the  kernel  is  a  nonempty  set  of 
spatial  discrete  locations  C  for  which  filter  values  can 
be  nonzero.  Except  for  locations  in  the  support  set, 
the  filter  value  is  0: 

fc(x,  y )  =  0,  (x,  y)  £  C.  (18) 

The  larger  the  filter  support,  the  better  the  perfor¬ 
mance,  but  small  kernels  can  be  highly  effective.26 

We  derive  the  optimal,  spatially  constrained  filter 
by  minimizing  MSE  €2  with  respect  to  the  elements  in 

C,  which  mathematically  requires  that 

mkrT0'  VM6C-  (19) 

These  constraints  can  be  expressed  in  a  system  of 
linear  equations26: 

2  d>p(x-x',  y-y%(x',  y')  =  ®sjr(x,  y), 

bc',y')EC 

V(x,  y)  G  C,  (20) 

where  'fy  is  the  autocorrelation  of  the  displayed  im¬ 
age  and  'fy.  p  is  the  cross  correlation  of  the  scene  and 
the  displayed  image.  The  number  of  equations  and 
the  number  of  unknowns  are  both  equal  to  the  num¬ 
ber  of  elements  in  support  set  C;  i.e.,  there  are  |  C  j 
equations  in  |  C  \  unknowns. 

C.  Parametric  Cubic  Convolution 

Piecewise  cubic  convolution  is  a  popular  interpola¬ 
tion  method  for  image  reconstruction  that  is  tradi¬ 
tionally  implemented  by  separable  convolution  with 
a  small  one-dimensional  kernel  consisting  of  piece- 
wise  cubic  polynomials.27-28  One  can  generalize  this 
popular  method  to  two  dimensions  and  reformulate 
it  by  relaxing  constraints  to  perform  reconstruction 
and  restoration  in  one  pass  with  small-kernel  con¬ 
volution.29  With  constraints  for  symmetry,  continu¬ 
ity,  and  smoothness,  the  two-dimensional  kernel 
with  support  [-2,  2]  X  [-2,  2]  has  five  parameters 
{oi,  a2,  a3,  a4,  a5}: 

fP(x,  y)  =fo(x,  y)  +  ai  /i(x,  y)  +  a2  f2(x,  y) 

+  a3  fa(x,  y)  +  a4  f4(x,  y)  +  a5  f5(x,  y), 

(21) 
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(g)  SA+Wiener 


(h)  Norm  2  Data 
Fig.  3.  Simulation  results. 


where  f0—f5  are  defined  in  Appendix  A.  This  kernel, 
designated  2D-5PCC-R  (to  distinguish  it  from  two- 
dimensional  piecewise  cubic  interpolation30),  is  a  con¬ 
tinuous  function. 

One  can  derive  the  optimal  2D-5PCC-R  kernel  f„ 
for  an  ensemble  of  scenes  by  minimizing  the  MSE  e2 
with  respect  to  the  five  parameters.  Computing  the 
partial  derivatives  of  e2  with  respect  to  the  parame¬ 
ters  and  solving  for  simultaneous  equality  with  zero: 


yields  five  equations  for  the  optimal  parameter  value: 

J  J  fi(u,  ^){Re[d>s,p(w,  iO]_/o(“>  v)%(u>  cOjdudu 

=  fi{u,  V)\fp{u,  v)-fo(u,  0)]4>p(a,  u)dwdu, 


de2 

de2 

de2 

de2 

d€‘ 

da1 

II 

1 

8 

ro 

II 

da3 

da4 

da. 

(22) 


i  =  1.  .  .5. 


(23) 
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(a)  Small-kernel  Wiener  filter 


(b)  2D-5PCC-R 

Fig.  4.  Small  reconstruction  and  restoration  kernels  for  the  sim¬ 
ulation  experiment. 

4.  Experimental  Results 

In  this  section  we  present  experimental  results  for  a 
simulated  imaging  system  and  for  real  images.  In  the 
simulation  the  scene  is  a  high-resolution  digital  im¬ 
age.  The  simulated  scene  is  blurred,  sampled,  and 
degraded  by  noise  (by  digital  processing)  to  produce 
simulated  microscanned  images.  These  images  are 
reconstructed  and  restored  by  the  optimal  CDC 
Wiener  filter,  by  the  small-kernel  Wiener  filter,  by 
the  2D-5PCC-R  filter,  by  shift-and-add  with  Wiener 
deconvolution,  denoted  SA  +  Wiener,20-31  and  by 
norm  2  data  with  L1  regularization  (denoted  Norm  2 
Data).23-31 

We  implement  the  superresolution  computations 
for  the  optimal  Wiener  filter  in  the  Fourier  frequency 
domain  by  multiplying  the  filter  defined  in  Eq.  (17)  by 
the  Fourier  transform  of  the  registered,  combined 


Fig.  5.  Low-resolution  infrared  image  of  a  four-bar  target  used  for 
estimating  the  acquisition  transfer  function. 

images  defined  in  Eq.  (5).  The  computations  for  the 
small-kernel  Wiener  filter  and  the  2D-5PCC-R  filter 
are  implemented  in  the  spatial  domain,  for  computa¬ 
tional  efficiency,  by  convolution  of  the  registered  com¬ 
bined  images  with  the  small  kernels  defined  by  the 
solutions  for  Eqs.  (20)  and  (23),  respectively.  Both  SA 
+  Wiener  and  Norm  2  Data  were  developed  at  the 
multidimensional  Signal  Processing  Research  Group 
at  the  University  of  California,  Santa  Cruz.31  Norm  2 
Data  is  an  iterative  superresolution  method,  and 
these  experiments  use  default  parameter  values  (ex¬ 
cept  the  deconvolution  kernel)  of  the  software. 

For  simulated  imaging,  all  resultant  images  are 
compared  to  the  original  scene.  Because  the  phase 
shifts  between  the  simulated  scene  and  the  mi¬ 
croscanned  images  are  known,  the  simulation  facili¬ 
tates  true  quantitative  measures  of  reconstruction 
and  restoration  performance.  We  also  apply  the  su¬ 
perresolution  methods  to  real  images  acquired  by 
panning  an  infrared  camera  slowly  across  a  fixed 
scene.  Because  the  true  scene  values  are  unknown, 
quantitative  evaluation  is  not  possible.  Also,  the  mi¬ 
croscanning  shifts  must  be  estimated,  so  the  results 
are  affected  by  registration  error.  Nonetheless,  real 
images  are  useful  for  qualitatively  demonstrating  the 
effectiveness  of  the  small-kernel  methods  in  practice. 

A.  Simulation  Results 

Figure  3(a)  illustrates  a  256  X  256  digital  image 
acquired  by  aerial  photography32  that  is  used  as  a 
simulated  scene,  and  it  is  therefore  the  ideal  super¬ 
resolution  image.  The  simulated  scene  is  blurred  by  a 
Gaussian  low-pass  filter  to  simulate  acquisition  blur¬ 
ring: 

h(u,  v )  =  exp[-(ii2  +  a2)],  (24) 

so  the  system  transfer  function  at  the  Nyquist  limit 
along  each  axis  is  /i(0.0,  0.5)  =  h{ 0.5,  0.0)  =  0.779. 


Table  1.  Fidelity  and  Computational  Costs  of  Various  Methods 


Performance  Metric 

Computation  Method 

CDC  Wiener 

Small-Kernel  Wiener 

2D-5PCC-R 

o-Moms 

SA  +  Wiener 

Norm  2  Data 

Fidelity 

0.980 

0.975 

0.966 

0.928 

0.966 

0.979 

Preprocessing  time  (s) 

0.521 

0.781 

0.940 

0.000 

0.150 

0.000 

Filtering  time  (s) 

0.121 

0.030 

0.030 

0.120 

0.380 

49.121 
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After  blurring,  sixteen  64  X  64  low-resolution  images 
are  created  with  simulated  microscanning  at  quarter- 
pixel  intervals  along  each  axis.  In  each  simulated 
microscanned  image,  Gaussian  white  noise  is  added 
to  each  pixel  such  that  the  blurred  signal-to-noise 
ratio  (BSNR)  is  30  dB: 

BSNR  =  10  logl0(o"p  2/o"e  2),  (25) 

where  up2  is  the  variance  of  the  blurred  microscanned 
images  (after  blurring  and  before  additive  noise)  and 
ue2  is  the  variance  of  the  noise.  Figure  3(b)  illustrates 
one  of  the  microscanned  images  (reference  image  p0) 
interpolated  back  to  256  X  256  resolution  by  nearest- 
neighbor  interpolation  to  show  the  granularity  of  the 
sampling.  Figure  3(c)  illustrates  a  higher-quality  in¬ 
terpolation  obtained  with  cubic  o-Moms.33  Four  of  the 
sixteen  microscanned  images  are  used  for  deriving 
the  reconstruction  and  restoration  filters  and  for  su¬ 
perresolution  processing: 


X 

0 

X 

0 

0 

0 

0 

0 

X 

0 

X 

0 

0 

0 

0 

0 

where  X  and  O,  respectively,  stand  for  locations  (at 
quarter-pixel  intervals)  with  and  without  samples. 

The  actual  scene  power  spectrum  is  used  to  derive 
the  optimal  CDC  Wiener  filter  to  benchmark  the  op¬ 
timal  fidelity.  The  power  spectrum  for  optimizing  the 
small-kernel  Wiener  filter  and  the  2D-5PCC-R  filter, 
however,  is  estimated  (as  typically  is  required  in 
practice)  using  the  power-spectrum  model  of  a  two- 
dimensional  isotropic  Markov  random  field  (MRF).34 
The  model  can  be  fitted  to  the  image  power  spec¬ 
trum30  or  interactively  parameterized  for  visual  qual¬ 
ity  (as  is  done  here).  The  MRF  autocorrelation  is 


<J>s(x,  y)  =exp(  — v'x2  +  y2/ p),  (26) 

where  p  is  the  mean  spatial  detail  of  the  scene  in  pixel 
units.  The  mean  spatial  detail  can  be  interpreted  as 
the  average  size  of  the  details  in  the  scene.  In  terms 
of  the  Hankel  transform,35  the  power  spectrum  of  the 
isotropic  MRF  is 


<Fs(w,  v) 


2Trp2 

[1  +  4tt2P2(w2  +  o2)]3/2- 


(27) 


The  CDC  Wiener  filter,  the  small-kernel  Wiener  filter 
with  support  limited  to  [-2,  2]  X  [-2,  2],  and  the 
2D-5PCC-R  filter  were  derived  for  this  simulation 
based  on  the  isotropic  MRF  scene  model  with  a  mean 
spatial  detail  of  4  pixels.  Figure  5  illustrates  the 
small-kernel  Wiener  PSF  [Fig.  4(a)]  and  the  2D- 


Fig.  6.  Superresolution  average  scan  of  the  bar  target. 


5PCC-R  PSF  [Fig.  4(b)].  The  optimal  parameter  val¬ 
ues  for  2D-5PCC-R  are  a:  =  74.176,  a2  =  -95.360, 
a3  =  16.804,  a4  =  —0.967,  and  a5  =  0.238. 

Figure  3(d)— 3(h)  presents  the  simulation  results  for 
the  superresolution  methods.  To  limit  boundary  ef¬ 
fects,  the  borders  of  all  resultant  images  are  cleared. 
Visually,  the  superresolution  images  produced  by  all 
the  restoration  and  reconstruction  filters  [Figs.  3(d) — 
3(h)]  are  better  than  provided  by  a  single  frame 
[Figs.  3(b)  and  3(c)].  For  example,  the  small,  lightly 
colored  rectangle  in  the  upper-left  quadrant  between 
the  diagonal  runway  and  the  leftmost  vertical  run¬ 
way  are  clearer  in  the  superresolution  images.  The 
images  produced  by  the  five  superresolution  methods 
are  of  similar  visual  quality,  but  the  CDC  Wiener 
filter  appears  to  produce  the  best  image  and  the  SA  + 
Wiener  filter  appears  to  produce  the  worst  image. 
The  image  produced  by  the  small-kernel  Wiener  filter 
appears  to  be  slightly  sharper  than  the  image  pro¬ 
duced  by  the  2D-5PCC-R  filter. 

Table  1  lists  the  quantitative  fidelity  and  compu¬ 
tational  costs  for  the  various  methods  with  the  sim¬ 
ulation.  The  computational  costs  were  measured  in 
seconds,  averaged  over  multiple  runs  by  MATLAB 


Fig.  7.  Estimated  acquisition  transfer  function  hx(u). 
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(a)  256  x  256  frame 


(b)  Cubic  o-Moms 


(c)  Small-kernel  Wiener  filter 


(e)  SA+Wiener  (f)  Norm  2  Data 

Fig.  8.  Superresolution  results  for  a  microscanned  infrared  system. 


6.5  Release  13  program  on  an  IBM  R32  (Intel  Pen¬ 
tium  M  1.8  GHz  CPU,  256  MB  RAM,  MS  Windows 
XP  Professional  2002).  The  optimal  CDC  Wiener  fil¬ 
ter  (which  uses  the  actual  scene  power  spectrum)  has 
the  highest  fidelity,  as  expected  mathematically.  The 


CDC  Wiener  filter  requires  preprocessing  for  comput¬ 
ing  the  filter  (which  then  can  be  used  to  filter  multiple 
images)  and  requires  forward  and  inverse  Fourier 
transforms  for  application  of  the  filter.  As  expected, 
the  image  from  the  cubic  o-Moms  with  a  single  frame 
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has  the  lowest  fidelity.  The  iterative  method,  Norm  2 
Data,  has  a  fidelity  nearly  equal  to  that  of  the  CDC 
Wiener  filter  but  requires  50  iterations  and  nearly 
50  s  for  Fig.  3(h).  Of  the  three  small-kernel  filters  that 
can  be  efficiently  applied  by  spatial  convolution,  the 
small-kernel  Wiener  filter  has  somewhat  higher  fi¬ 
delity  than  the  2D-5PCC-R  or  the  SA  +  Wiener  filter. 

The  visual  and  quantitative  results  for  the  simu¬ 
lation  indicate  that  the  small-kernel  Wiener  filter 
and  2D-5PCC-R  effectively  improve  image  quality 
with  efficient  spatial-domain  processing. 

B.  Results  for  Real  Images 

The  superresolution  reconstruction  and  restoration 
methods  require  characterizations  of  the  system,  in¬ 
cluding  the  noise  power  spectrum  and  acquisition 
transfer  function.  Noise  can  be  characterized  accu¬ 
rately  from  flat-field  calibration  images  or  from  re¬ 
gions  of  uniform  background  of  acquired  images.  The 
transfer  function  can  be  estimated  to  frequencies  be¬ 
yond  the  Nyquist  limit  by  use  of  a  knife-edge  tech¬ 
nique  with  images  for  various  sample-scene  phase 
shifts.36  Microscanning  provides  such  images. 

For  this  experiment  with  real  images,  the  acquisi¬ 
tion  transfer  function  of  an  infrared  camera  system 
was  estimated  from  microscanned  images  of  a  four-bar 
target.  The  camera  platform  was  microscanned  as  low- 
resolution  images  were  recorded.  Figure  5  illustrates  a 
small  piece  of  one  of  the  256  X  256  low-resolution 
images  from  the  microscanned  sequence.  Employing 
the  superresolution  knife-edge  technique,36  we  ob¬ 
tained  and  averaged  a  sequence  of  120  low-resolution 
images  of  the  bar  target  with  subpixel  accuracy  (to  0.25 
pixel).  Figure  6  illustrates  the  resultant  superresolu¬ 
tion  horizontal  slice  across  the  four-bar  targets,  super¬ 
imposed  upon  the  model  of  the  bar  target  scene 
estimated  by  thresholding  of  the  registered  slice. 

Figure  7  illustrates  the  modulation  transfer  func¬ 
tion  estimated  by  the  CDC  Wiener  filter.  For  this 
experiment  the  two-dimensional  acquisition  transfer 
function  was  modeled  as  the  separable  product  of  the 
one-dimensional  estimate: 

h(u,  v)  =  hx(u)hx(v).  (28) 

A  more  accurate  estimate  of  the  two-dimensional  ac¬ 
quisition  transfer  function  could  be  made  from  two  or 
more  slices. 

In  this  experiment,  the  goal  of  superresolution  recon¬ 
struction  and  restoration  is  to  generate  a  1024  X  1024, 
high-fidelity  image  from  multiple  microscanned  256 
X  256  images.  Figure  8(a)  illustrates  one  of  a  mi¬ 
croscanned  sequence  of  low-resolution  images  from 
the  infrared  camera  system,  interpolated  to  1024 
X  1024  by  nearest-neighbor  interpolation.  The  scene 
is  modeled  as  a  modulation  transfer  function  with  a 
mean  spatial  detail  equal  to  8  pixels.  The  blurred 
signal-to-noise  ratio  was  estimated  to  be  30  dB. 

Figure  9  illustrates  the  small-kernel  Wiener  and 
cubic  convolution  reconstruction  and  restoration  ker¬ 
nels,  based  on  these  system  characterizations.  The  op¬ 
timal  parameter  values  for  2D-5PCC-R  are  oq  = 


(b)  2D-5PCC-R 

Fig.  9.  Small  reconstruction  and  restoration  kernels  for  the  real 
image  experiment. 


—  1.052,  a2  =  —7.033,  a3  =  7.584,  a4  =  —1.123,  and 
as  =  0.137.  The  small-kernel  Wiener  PSF  provides 
more  sharpening. 

Figure  8(b)  illustrates  a  1024  X  1024  image  recon¬ 
structed  from  a  single  image  by  cubic  o-Moms. 
Figures  8(c)-  8(f)  illustrate  1024  X  1024  images  recon¬ 
structed  and  restored  with  the  small-kernel  Wiener 
filter,  2D-5PCC-R,  SA  +  Wiener,  and  Norm  2  Data. 
These  results  were  generated  from  just  three  256 
X  256,  microscanned  low-resolution  images  with  dif¬ 
ferent  relative  shifts.  The  relative  shift  pattern  for  the 
three  images  is 


X 

X 

0 

X 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

To  limit  boundary  effects,  the  borders  of  the  pro¬ 
cessed  images  are  cleared.  Visually  the  resultant  im¬ 
ages  in  Figs.  8(c),  8(d),  and  8(f)  from  superresolution 
reconstruction  and  restoration  are  substantially  better 
than  a  single  image  in  either  Fig.  8(a)  or  Fig.  8(b). 
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Figure  9(e)  from  the  SA  +  Wiener  superresolution 
method  is  blurred.  The  superresolution  image  from  the 
small-kernel  Wiener  filter  is  slightly  sharper  than  the 
image  from  2D-5PCC-R,  consistent  with  the  simula¬ 
tion  results. 

5.  Conclusions 

Images  with  superresolution  and  high  fidelity  are 
required  in  many  imaging  applications.  In  this  paper 
we  have  investigated  superresolution  methods  for 
microscanning  imaging  systems  that  efficiently  im¬ 
plement  reconstruction  and  restoration  with  small 
convolution  kernels.  The  approach  is  based  on  an 
end-to-end  system  model  that  accounts  for  the  fun¬ 
damental  trade-off  in  imaging  system  design  between 
acquisition  and  aliasing.  The  system  model  also  ac¬ 
counts  for  noise,  misregistration  of  the  low-resolution 
images,  and  the  display  system.  Computationally  ef¬ 
ficient  reconstruction  and  restoration  are  achieved  by 
use  of  one-pass  convolution  with  small  kernels. 

We  have  developed  two  small  convolution  kernels 
for  improved  resolution  and  fidelity:  the  spatially  con¬ 
strained  Wiener  filter  and  a  parametric  cubic  convolu¬ 
tion  (designated  2D-5PCC-R).  Subject  to  constraints, 
both  have  been  optimized  with  respect  to  maximum 


end-to-end  system  fidelity.  Experimental  results  for  a 
simulated  imaging  system  and  for  real  images  indicate 
the  effectiveness  of  the  small-kernel  methods  for  in¬ 
creasing  resolution  and  fidelity.  Visually,  the  super¬ 
resolution  images  from  the  small-kernel  Wiener  filter 
are  slightly  sharper  than  images  from  2D-5PCC-R,  but 
both  small-kernel  methods  yield  significant  quantita¬ 
tive  and  qualitative  improvements. 

Additional  work  to  develop  efficient  implementa¬ 
tions  of  the  small-kernel  superresolution  methods  is 
required.  Even  with  the  efficiency  of  one-pass  resto¬ 
ration  and  reconstruction  by  use  of  small  kernels, 
superresolution  requires  significant  processing.  Each 
low-resolution  image  pixel  in  the  kernel’s  region  of 
support  around  each  superresolution  pixel  contrib¬ 
utes  to  the  output  value.  Superresolution  processing 
of  K  images,  with  kernel  support  of  [— S,  S]  X  [-S, 
S]  pixels  and  a  superresolution  increase  of  R  X  R, 
requires  4 KS2R2  floating-point  multiplications  and 
additions.  If  the  distribution  of  low-resolution  pixels 
varies  with  respect  to  the  superresolution  pixels, 
multiple  kernels  (each  with  KS2  weights)  should  be 
used  in  the  computation.  These  issues  make  imple¬ 
mentation,  especially  in  hardware,  challenging. 


Appendix  A.  Components  of  the  2D-5PCC-R  Cubic  Convolution  Kernel  for  Superresolution.  Restoration  and 
Reconstruction  [Eq.  (21)] 


fo(x,  y )  = 


fi(x,  y )  = 


f2(x,  y)  = 


fa(x,  y)  = 


fi(x,  y)  = 


fs(x,  y)  = 


r  2  2  2  2  i  i 

x y  -  x  —  y  +1 

0  <x  <  1, 

0  <y  <  1, 

(2xy2  —  2x  —  2 y2  +  2)(x  —  2)2 

1  <  x  ^  2, 

0  <y  <  1, 

(4xy  —  4y  —  4x  +  4)(x  —  2)2(y  —  2)2 

1  <  x  £  2, 

1  <y  <2, 

^(2x2y  -2 y  -  2x2  +  2 )(y  -  2)2 

0  <  x  <  1, 

l<y  <2, 

'x3y3  -  x2y2 

0  <  x  <  1, 

0  <y  <  1, 

(5xy3  -  4 xy2  -  4 y3  +  3 y2)(x  -  2)2 

1  <  x  <  2, 

0  <y  <  1, 

(9xy  —  8y  —  8x  +  7)(x  —  2)2(y  —  2)2 

1  <  x  <  2, 

l<y  <2, 

^(5x3y  —  4x2y  —  4x3  +  3  x2)(y  —  2) 2 

0  <  x  <  1, 

l<y  <2, 

x3y2  -  2x2y2  +  x2y3 

0  <  x  <  1, 

0  <y  <  1, 

(4xy3  -  3xy2  -  3y3  +  2 y2)(x  -  2)2 

1  <  x  <  2, 

0  <y  <  1, 

(8xy  —  7y  —  7x  +  6)(x  —  2)2(y  —  2)2 

1  <  x  <  2, 

1  <y  <2, 

^(4x3y  -  3 x2y  -  3x3  +  2x2)(y  -  2) 2 

0  <  x  <  1, 

1  <y  <  2, 

"x3  -  x2  -  y2  +  y3 

0  <  x  <  1, 

0  <y  <  1, 

(2xy3  -  2xy2  +  x  +  y2  -  y3  ~  l)(x  -  2)2 

1  <  x  <  2, 

0  <y  <  1, 

(4xy  —  3x  —  3y  +  2)(x  —  2)2(y  —  2)2 

1  <  x  <  2, 

l<y  <2, 

^(2x3y  -  2x2y  +  y  +  x2  -  x3  -  l)(y  -  2)2 

0  <  x  <  1, 

1  <y  <2, 

r5x2  -  6x2y2  +  5 y2  -  4 

0  <  x  <  1, 

0  <y  <  1, 

(l3y2  —  14xy2  +  12x  —  ll)(x  —  2)2 

1  <  x  <  2, 

0  <y  <  1, 

(30x  +  30y  —  32xy  —  28)(x  —  2)2(y  —  2)2 

1  <  x  <  2, 

l<y  <2, 

(13x2  -  14x2y  +  12 y  -  11)  (y  -  2)2 

0  <  x  <  1, 

1  <y  <  2, 

'Ax2  -  3 x2y2  +  4y2  -  4 

0  <  x  <  1, 

0  <y  <  1, 

(5 y2  —  Axy2  +  8x  —  8)(x  —  2)2 

1  <  x  <  2, 

0  <y  <  1, 
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1  <y  <2, 
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